using System;
using System.Reflection;
using NBody;
using NBody.Analysis;
using OptionParser;

[assembly: AssemblyTitle ("NBodyAngularMomentum")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Compute various angular momentum quantities.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class JOptions : InputOutputGetOptions
{
    [Option ("-b", "Number of bins")]
    public int NumBins;
    
    [Option ("-s", "Bin spacing")]
    public Bin.Spacing Spacing;

    [Option ("-l", "Log plot")]
    public bool LogPlot;
    
    [Option ("-sl", "Semi Log plot")]
    public bool SemiLogPlot;
    
    [Option ("-c", "Adjust for centre of mass")]
    public bool CentreOfMass;
    
    public JOptions()
    {
        NumBins = 100;
        Spacing = Bin.Spacing.EqualMass;
    }
}

public class NBodyAngularMomentum
{
    public static void Main(string[] args)
    {
        JOptions opts = new JOptions();
        opts.ProcessArgs(args, ".jbin");
                        
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        double jx, jy, jz;
        CalculateTotalJVector(s, out jx, out jy, out jz);
        s.AlignWithVector(jx, jy, jz);
        
        if (opts.CentreOfMass)
            s.AdjustForCentreOfMass(1e-3);
        
        Bin jbin = new AngularMomentumBin();
        jbin.BinSpacing = opts.Spacing;
        jbin.NumBins = opts.NumBins;
        jbin.Analyse(s);
        
        if (opts.LogPlot)
        {
            jbin.Data.Print(opts.OutFile, new DataFile.Transform(Math.Log10));
        }
        else if (opts.SemiLogPlot)
        {
            DataFile.Transform[] trans = new DataFile.Transform[jbin.Data.Columns];
            trans[0] = new DataFile.Transform(Math.Log10);
            for (int i = 1; i < jbin.Data.Columns; i++)
                trans[i] = new DataFile.Transform(DataFile.Table.IdentityTransform);
            jbin.Data.Print(opts.OutFile, trans );
        }
        else
            jbin.Data.Print(opts.OutFile);
    }
    
    static void CalculateTotalJVector(NBodySystem s, out double jx, out double jy, out double jz)
    {
        jx = jy = jz = 0.0;
        foreach (Particle p in s)
        {
            jx += p.Mass * (p.Y * p.Vz - p.Z * p.Vy);
            jy += p.Mass * (p.Z * p.Vx - p.X * p.Vz);
            jz += p.Mass * (p.X * p.Vy - p.Y * p.Vx);
        }
    }
}
